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ABSTRACT 

Surprisingly little is known about the origin and evolution of the Milky Way's satellite 
galaxy companions. UV photoionisation, supernova feedback and interactions with the 
larger host halo are all thought to play a role in shaping the population of satellites 
that we observe today, but there is still no consensus as to which of these effects, 
if any, dominates. In this paper, we revisit the issue by re-simulating a Milky Way 
class dark matter (DM) halo with unprecedented resolution. Our set of cosmological 
hydrodynamic Adaptive Mesh Refinement (AMR) simulations, called the Nut suite, 
allows us to investigate the effect of supernova feedback and UV photoionisation at 
high redshift with sub-parsec resolution. We subsequently follow the effect of inter- 
actions with the Milky Way-like halo using a lower spatial resolution (50pc) version 
of the simulation down to z = 0. This latter produces a population of simulated 
satellites that we compare to the observed satellites of the Milky Way and M31. We 
find that supernova feedback reduces star formation in the least massive satellites but 
enhances it in the more massive ones. Photoionisation appears to play a very minor 
role in suppressing star and galaxy formation in all progenitors of satellite halos. By 
far the largest effect on the satellite population is found to be the mass of the host 
and whether gas cooling is included in the simulation or not. Indeed, inclusion of gas 
cooling dramatically reduces the number of satellites captured at high redshift which 
survive down to z = 0. 



1 INTRODUCTION 



' In the standard cold dark matter paradigm of galaxy for- 
mation, galaxies grow inside dark matter halos that merge 
hierarchically. In other words, smaller halos are captured by 
larger halos and the galaxies they contain become satellite 
galaxies of the host galaxy until dynamical friction finally 
forces them to coalesce. Early attempts to reproduce the 
observed Milky Way satellite population using dark-matter- 
only simulations overproduced the number of low mass satel- 
lites by several orders of magnitude when attributing to 
each simulated dark mat ter satellite a gala xy using a con- 
stant mass-to- light ratio (|Moore et al.lll99 9l). Current state- 
of-the-art dark mat ter simulations, such a s the Aquarius 
(ISpringel et aljfeOQSl ) and the Via Lactea II (|Diemand et al.l 
l2008h still reach the same conclusion. 

Whilst the existing sample of Mil ky Way satellite galax- 
ies is almost certai nly incomplete (|Koposov et al.l IboOSt 
iTollerud et ail jflflj ) and new satellite galaxie s are con- 
tinually being discovered (e.g. iBelokurov et al.l (|2010h h it 
appears extremely unlikely that new observations will un- 
cover galaxies populating every dark matter substructure 
predicted to exist around the Milky Way. As a result, vari- 
ous authors have attempted to explain this discrepancy by 
invoking physical mechanisms that reduce or prevent star 



formation in the majority of the smaller halos, making these 
dark matter substructures fainter or completely dark in the 
process. In particular, photoionisation and feedback from su- 
pernovae have been proposed as the most likely mechanisms 
to prevent gas from condensing to form stars, although other 
mechanisms, such as cosmic rays, have also been suggested 
jWadepuhl fe Springelllioiol ). 

The UV ionising background has been argued to 
be effective at halting or preventing star formation in 
low-mass halos in studies using analytic arguments, ob- 
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Gnedinl (|2000T ), hydrodynamic simulations of Milky Way- like 
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dHoeft et al.l 120061; iLibeskind etail 120101; iNickerson et al.l 
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Many of these authors have found a best fit to the ob- 
served luminosity function by adopting instantaneous 
reionisation at z ~ 11, in agreement with the latest 
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WMAP e stimate (z re i on = 1 5 ± 1.2 jLarson et alJkoid)) . 
However, IHoeft et all (|2006h : IWadepuhl fc Springell < |201Ch ; 
iGuo et alJ (|2010l ) also find reionisation to have limited 
effectiveness in completely suppressing star formation in 
low-mass h a los t hat have already begun forming stars. 
IHoeft et all |2006h determine that there is a characteristic 
mass of 6.5 x 10 9 /i _1 Mq below which halos become unable 
to ret ain baryons down to z = 0, and hence cannot form 
stars. lOkamoto et al.l (|2008l ) find a similar result, which 
they translate into a minimum circula r velocity v m ax of 
25km/ s below which halos are dark. lOkamoto fc Frenkl 
|2009l ) further refine this finding, stating that the cut-off 
should be lowered to v ma x = 12km/s at reionisation (z = 9 
in their case). To further obscure the picture, the extent of 
the epoch of reionisation itself is poorly constrained, with 
only a lower limi t of z > 6 on its completion provided by 
observ ations (e.g. ICen et al.l (|2009l ); iMesinger fc Furlanettd 
(HI)). 

Another process is proposed in iDekel fc Silk! (l986). 
They argued that the suppression of star formation by 
supernova feedback in dwarf galaxies embedded in dark 
matter halos could explain the obser ved scaling relation s 
in luminosity , metallicity and radius |Dekel fc Wool 120031 ) . 
I Benson et al.l (|2003l ) suggested that supernovae could help 
explain the unexpectedly low dwarf galaxy luminosity func- 
tion. By removing the gas from galaxies via the injec- 
tion of thermal and kinetic energy into the interstellar 
medium, supernovae would reduce the n umber of stars 
forme d i nside dwarf halos. A uthor s s uch as Low fc Ferraral 
d 19991); [Mashchenko et all (|2008l ); iRicotti et all (|2008l ); 
ICeverino fc Klypinl l|2009l ) manage t o generate massive 
supernova-driven galactic winds, but iTassis et al.l (|2008l ) 
claim that including supernova feedback does not affect the 
properties of their simulated galaxies, attributing the scal- 
ing relations found in dwarf galaxies to low star formation 
efficiencies in weak potentials instead. These contradictions 
find their origin in the different numerical recipes and nu- 
merical resolutions adopted, along with the variety of galaxy 
masses and merger histories used. As a consequence, it is 
still unclear as to precisely what effect supernovae have on 
the ISM gas, and hence on the star formation in low-mass 
galaxies. Indeed, supernovae are potentially able to drive 
either positive or negative feedback cycles. Outflows from 
supernovae can remove gas from the galaxy, preventing it 
from forming stars. However, they also release metals into 
the surrounding ISM; metal line cooling increases the effi- 
ciency of gas cooling and hence p romotes the collaps e of gas 
clouds into star-forming regions (|Powell et al ] |201ll ). More- 
over, blast wave com pression has been obser ved to trigger 
star-forming regions (jAssousa fc Herbst|[l98(il ). 

It has also been recently (re-) suggested that ou r mode l 
of dark matter is incorrect. iBovlan-Kolchin et al.l (|201lT ); 
lLovell et all (|2012h ;l Di Cintio et al.l (|201lM examine the rela- 
tionship between the maximum circular velocity of the dwarf 
spheroidal satellites and the radius at which this velocity is 
found, and determine that the velocities found are higher 
than those found in either ACDM pure dark matter simula- 
tions or simulations with baryon physics using smoothed 
partic le hydrodynamics (SPH). In fact, |Pi Cintio et al.l 
|201ll ) determine that gas cooling makes the problem worse, 
since the central density of the halo increases and the radius 
at which the maximum circular velocity is found decreases. 



lLovell et all (|2012h suggests that warm dark matter (WDM) 
resolves the problem, but they do not run simulations with 
both WDM and baryon physics. 

Bearing these caveats in mind, the present paper aims 
to constrain satellite galaxy formation and evolution, and 
more specifically the role played by supernova feedback and 
reionisation in the process. To this end, we use the ^zz*& 
(nut) suite of high resolution hydrody namic cosmological 
simulations of a M ilky Way- like galaxy (|Powell et al.ll201ll ; 
iKimm et al.ll20l"ll ). At scales of l-10pc, these resolve large 
molecular clouds, and hence model the interstellar gas and 
stellar feedback in grea ter detail, which potentially affects 
star formation histories (Slvz et al 120051 ) . This approach dif- 
fers from previous work investigating the Milky Way satel- 
lites using hydrodynamic simulations as these generally cap- 
ture the ISM at lower resolution, and attempt to compensate 
for this by introducing analytic e xpressions to account for 
the multiphase ISM and o utflows (|Scannapieco et ai1l2005l . 
120061 ; iMurante et al.ll201ol) . As it is currently too costly to 
simulate a Milky Way-sized halo and its substructures with 
parsec resolution and hydrodynamics throughout the life- 
time of the Universe, most of our analysis is restricted to 
high redshift (z > 6). However, observational studies con- 
clude that a vast majority of satellite galaxies contain stars 
which formed prior to z ~ 2 an d in many ca s es pri or to 
even z ~ 5 (see r ecent review bv lTolstov etZI (11)091 ) and 
also iKirbv et al.l (|201lT )). Therefore, a high redshift study 
of these objects should be able to shed light on the prob- 
lem, provided one is able to accurately predict their spatial 
distribution at z = 0. 

This latter requirement is in itself a major challenge as 
it presumably requires hydrodynamics simulations which in- 
clude (at least) radiative cooling. Indeed, since gas cooling 
can significan tly increase the central density of dark mat- 
ter halos (e.g. iBlumenthal et al.l (|l984f )). one expects physi- 
cal processes like dynamical friction and tidal disruption of 
the satellites to be altered. The extent of these differences 
needs to be quantified because most of the studies men- 
tioned earlier in this introduction rely on pure dark matter 
simulations to underpin analytic arguments or graft semi- 
analytic models of galaxy formation. Several groups have 
looked at differences between simulations of galactic halos 
containing baryons and their pu r e dark matter N-body coun- 
terparts. For instance IPeiran 3 (|201Ch found that identical 
simulations of a local group-like volume with and without 
baryons matched well, but did not comment on satellites 
within halo s. In their cons t rained simulations of the Lo- 
cal Group, iLibeskind et al.l (|2010l ) find more satellite ha- 
los when baryons are included than in the identical pure 
dark matter run. Their radial distribution is also significantly 
more concentrated. By contrast, although they also follow a 
more concentrated radial distribution, satellites in t he bary- 
onic simulations of iRomano-Di'az et al.l (|2010l . 12009] ) survive 
for shorter times than their pure dark matter counterparts, 
which yields an overall lowe r number of satellites with i n R v j T 
when baryons are included. Schewtschenko fc Maccicl (|201ll ) 
finds a similar result to ILibeskind et al.1 (|2010l )~ in terms of 
number of halos but sugg ests that these results are n ot in- 
com patible with those of Romano-Di'az et al.l (|2010l ). Fur- 
ther, iD'Onghia et al.l (|201(]| ) find that the presence of a disk 
can affect the mass function of satellites around a host halo. 
We re-examine this issue using Eulerian AMR grid hydro- 
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dynamics instead of Lagrangian smoothed particle hydro- 
dynamics (SPH), and improving on both mass and force 
resolution for the dark matter. 

This paper is split into three main parts. In section[2] we 
describe the simulations, our algorithms for comparing them 
and for tracking halos at high redshift down to z — 0. Sec- 
tion [3] looks at the effect of feedback mechanisms on satellite 
galaxy formation at high redshift and ultra-high resolution. 
The third part, section SI is devoted to the present epoch, 
and how our results at high redshift affect the satellite pop- 
ulation we see today. 



2 METHODS 

In this section we discuss the methods employed to carry 
out the simulations used in this paper and the subsequent 
analysis techniques. 



2.1 Numerical simulations 

We analyse five simulations in the ^ = ^3 (nut) suite of sim- 
ulations (i=i=oiS is t he Ancient Egyptia n goddess of the sky) 
(|Powell et al.ll201ll ; iKimm et allfeullD . S = «sll is a set of 
cosmological hydrodynamic resimulations of a Milky Way- 
like halo at z — (throughout this paper, this halo will 
be referred to as the "Milky Way"). To run these simula- 
tions, we use the Adap tive Mesh Refinement (AMR) code 
RAMSES l|Tevssier]|2002l ). Each simulation starts from iden- 
tical initial conditions, which are gener ated with mpgrafic 
|Bertschingeril200ll ; iPrunet et aflfeoOSt ) using cosmological 
parameters consisten t with the WMAP 5 year measurements 
(|Dunklev et al.ll2009l ). The simulation volume is a periodic, 
cubic box of length 9h _1 Mpc with a minimum resolution 
of 128 3 dark matter particles and the same number of grid 
cells. Within this volume we carve out a spherical region of 
radius 1.44h _1 Mpc, centred on a halo that reaches a virial 
mass M v i r = 5 x 10 11 Mq at z=0. We place three nested 
grids in this spherical region with effective resolutions of 
256 3 , 512 3 and 1024 3 dark matter particles and grid cells. 
The minimum dark matter particle mass inside this region is 
equal to 5.6xlO 4 M0 (with the exception of the Dark Mat- 
ter run, which we describe later in this section). We then 
allow the grid inside the refinement region to adaptively re- 
fine up to a given maximum level for each simulation. Our 
refinement strategy is quasi-Lagrangian: a grid cell is refined 
when the number of dark matter particles in the cell exceeds 
8, or the baryonic mass of the cell reaches 8 mspH, where 
mspH = 9.4 x 1O 3 M0. The simulation parameters used are 
summarised in table [I] and in the text below. 

The three main simulations that we consider in this 
paper contain dark matter, gas cooling and a uniform UV 
background switched on at z = 8.5 to model reionization 
lHaardt fc Madaul (| 19961 ). Gas cooling is modelled as radia- 
tive energy loss from atomic processes including emission 
line cooling (below 10 4 K), with a primordial metallicity of 
0.001 Zq. Star formation in the simulation proceeds ac- 
cording to a Schmidt law on a local dynamical timescale 
|Cen fc Ostrikerll 19921 ) with an efficiency of 1%. The density 
threshold for star formation is set in each simulation to be 
comparable to the corresponding Jeans density pj of a cell 
on the highest level of refinement with a temperature of 100 



K. pj is given by (7rcf) / '(AjG) = ksT / { mu AjG) for an ideal 
gas, where Aj is set to the cell length (|Binnev fc Tremainel 
l2008h . 



We first 
erence run" . 
adaptively to 



run a simulation that we call the "Ref- 
The Reference run is allowed to refine 
up to 8 times inside the fixed refinement 
region, such that the densest regions are allowed to 
reach a maximum physical resolution of 50pc at all 
times, between a few times and an order of magnitude 
higher than other cosmological hydrodynamics simu- 
lations of Milky Way satellites ilLibeskind et al.l [20T0I; 
Nickerson et al.ll201ll; lOkamoto fc Frcnk 20 091; iParrv et al.1 



ISawala et al.l 



~l20Tj; 



Ricotti fc Gnedi nl |2005l : iRomano-Di'az et al 



Scannapicco ct al. 



Schewtschenko fc Macciol l201ll ; IWadepuhl fc Springell 
2010h . The minimum star particle mass in this simulation 



2010, 
201 ll; 



is 3.5xl0 4 M©. We run this simulation to z = 0. Note 
that the main purpose of this run is to act as a lower 
spatial resolution "twin" of the two "high resolution" 
simulations in this study, allowing us to determine which 
of the galaxies formed at high redshift are progenitors of 
Milky Way satellite galaxies today. For this reason the DM 
mass resolution is kept identical in all runs. The threshold 
for star formation in the Reference run is 10 atoms/cm 3 . 

We then run two high resolution simulations which are 
allowed to refine adaptively by up to 15 times so that its 
physical spatial resolution in the densest regions can reach 
a maximum of 0.5pc at all times. The first of these high 
resolution simulations we call the "Cooling run". As with 
the Reference run, the Cooling run contains dark matter, 
gas cooling and a uniform UV background switched on at 
z — 8.5, but now the threshold density for star forma- 
tion is 10 5 atoms/cm 3 . As a result, the minimum star par- 
ticle mass formed in the Cooling run is 167M0. The sec- 
ond of the high resolution simulations is called the "Feed- 
back" run. The Feedback run is identical to the Cooling run, 
except that it also includ es supernova feedback. Following 
iDubois fc Teyssierl (|2008l ). supernovae are implemented as 
Sedov blast waves with a radius of 2 grid cells (lpc) around 
a star particle 10 Myr after it formed. Note that while we do 
not resolve individual stars, at this mas s resolution an d as- 
suming a Salpeter initial mass function (|Salpeterlll955l ). we 
get one supernova per star particle. We assume supernova 
events entrain 50% of the in itial mass of the star par ticle 
(ffw = 1 in the notation of IDubois fc Teyssierl (|2008l )) in 
a wind and have a metal yield of 0.1. The energy released 
is given by ^sN^^-esN, where m* is the mass of the star 
particle, msN and esN are typical values for a massive star 
going supernova and tjsn is the fraction of the total mass 
formed that is turned into supernova ejecta. For this simu- 
lation, we use 7?sn = 0.106, msN = IOM and esN = 10 51 ergs 
ijPowell et al.ll201ll : lK"imm et al.ll201ll ). This translates into 
a minimum star particle final mass of 76M0 for the Feedback 
run. 

The density threshold for the Reference run is cho- 
sen to best match the star formation rate per halo mea- 
sured in the Cooling and Feedback runs, since these cap- 
ture the length scales of molecular clouds, allowing for a 
more realistic model of star formation. The density thresh- 
old has to obey two constraints: (i) the star formation den- 
sity threshold should be smaller than the corresponding 
/0j(~ 40 at/cm 3 ) on the highest level of refinement and (ii) 
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Simulation 




m DM 


Rmu (level) 


Gas Cooling 


m„ 


SNe 


UV 


Reference Run 





5.6xl0 4 M Q 


50pc (18) 


/ 


3.5xl0 4 M Q 




/ 


Cooling Run 


6.7 


5.6x10 4 M o 


0.5pc (25) 


/ 


167M Q 




/ 


Feedback Run 


6.7 


5.6xlO 4 M 


0.5pc (25) 


/ 


76M 


/ 


/ 


Dark Matter Run 





6.7xlO 4 M 


50pc (18) 










Adiabatic Run 





5.6x10 4 M o 


50pc (18) 








/ 



Table 1. Table of properties of numerical simulations included in this paper. The columns are, from left to right, the simulation name, 
the lowest redshift reached by the simulation, the minimum dark matter particle mass, the maximum spatial resolution of the AMR 
grid (with the associated level of refinement in brackets), whether the simulation includes gas cooling, the minimum star particle mass, 
whether the simulation includes supernova feedback, and whether the simulation includes a UV background. The bottom two simulations 
are considered only in section For a complete description of the simulations, see section [2. II 



stars should not form in smooth filaments, which yields a 
lower bound on the density t hreshold that we e mpirically 
determine to be ~ 10 at /cm 3 jPowell et al.ll201ll ). A higher 
threshold will thus limit star formation to cells with a higher 
average density. Note, however, that this is the average den- 
sity, and a volume of the ISM of length 50pc with a low 
average density may still host small regions of high density 
gas. Hence a gas cell in the Reference Run with a density 
below pj may still form stars. We determine that a value of 
10 atom/cm 3 better matches the star formation rates found 
in the Cooling and Feedback runs. This value was found 
by using a number of test runs of the Reference Run using 
different star formation density thresholds. 

Finally, we perform a further two simulations. These 
are called the Dark Matter run and the Adiabatic run. Both 
have the same initial conditions and refinement criteria as 
the Reference run. The Dark Matter run is a pure N-body 
dark matter simulation, in which the mass in baryons is 
replaced by mass in dark matter, such that a dark matter 
particle is 1/(1 — ft) times the mass of a dark matter particle 
in the runs containing baryons, where i s the universal 
baryo n fraction (0.173, based on the data in iDunklev et all 
(200!))). This gives it a minimum dark matter particle mass 
of 6.7xl0 4 M o rather than 5.6xlO 4 M , which is the value 
common to all the other runs. The Adiabatic run is identical 
to the Reference run, except that the gas is not permitted to 
radiate away its energy. As a result, no star formation takes 
place in the Adiabatic run, though we still include the UV 
background for sake of comparison. These two simulations 
are used to determine the effect of including more physics 
on satellite galaxy evolution from z ~ 6 (the redshift where 
the high resolution simulations stop) to z = 0. We discuss 
the results of this study in section POl and compare the Dark 
Matter run to other pure N-body dark matter simulations 
of Milky Way satellites in section 14.21 



2.2 Halo identification 

We use HaloMaker to identify dark matter halos and galaxies 
in each simulati on output using th e Most Massive Subhalo 
Method (MSM) (|Tweed et alj|2009h . We define independent 
halos as dark matter overdensities not contained within an- 
other halo, and subhalos as halos that are identified as sub- 
structures of other halos. Similarly, we define galaxies as 
overdensities in the star particles. An overdensity is defined 
as a structure which is above 178p cr it, where p cr u is the crit- 
ical density of the universe. The method works as follows. 
The density of each particle is found using the SPH tech- 



nique (e.g. lSpringel et all |200ll )). Peaks in the density field 
are then identified, which correspond to the centre of a halo 
or galaxy. Particles are then attached to a given peak by 
stepping through decreasing density thresholds, and assign- 
ing each particle above this threshold to the nearest peak. 
Saddle points in the density field are identified, which are 
used to construct a tree of peaks, truncated at 178p cr it- Each 
leaf of this tree is a halo or subhalo. The host halo is identi- 
fied as the most massive peak, while the other peaks become 
subhalos. For this procedure, we reject any identified halo 
that contains less than 40 particles, twice the absolute min- 
im um threshold befo re spurious halos are identified given 
in (|Tweed et alj|2009l )). We also reject any galaxy that con- 
tains less than 10 particles - this figure is lower because stars 
are typically more tightly clustered than dark matter halos). 
The minimum total mass of a given dark matter (sub) halo is 
thus 2.2xl0 6 MQand the minimum stellar mass for a galaxy 
is 3.5xlO 5 M in the Reference run, 17OOM in the Cooling 
run and 76OM in the Feedback run. 

We perform this process in every simulation for both 
the dark matter and the stars whenever possible. We thus 
identify every dark matter halo and every galaxy above the 
mass limits given in the last paragraph. We visually inspect 
the results of the halo identification by overplotting the halos 
on a projection of the density field and tune the halofinder 
parameters such that we minimise spurious halo detection 
or halos that are visually identifiable but not detected by 
the halo finder. 

We then identify the host halo of galaxies to determine 
which halos are luminous and which are dark. This sorting 
proceeds in two steps: (i) for each galaxy, we make a list of 
halos that enclose it within their virial radius (ii) we select 
the halo that lies closest to the galaxy centre. This final step 
is needed when the galaxy's host is a subhalo of a larger 
halo, and the galaxy lies within the virial radius of both the 
subhalo and its host. This allows us to match halos between 
runs and to compare the properties of the embedded galaxies 
in each simulation on an individual halo basis. 



2.3 Halo twinning 

In order to compare the simulations, we associate each halo 
in the Reference run with a counterpart (called 'twin') in the 
Cooling and Feedback runs. As well as allowing us to com- 
pare halos between different runs at common redshifts, this 
procedure also permits us to track the halos in the Cooling 
and Feedback runs down to z = via the merger trees of 
their counterparts in the Reference run (see section [2~4)l . 
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Figure 1. Visualisation of the halo twinning results at z = 6.7, showing all successfully twinned halos. The left panel shows the cubic 
volume (300 physical kpc on a side) containing all the Milky Way satellite progenitor dark matter halos identified by z = 0. The right 
panel shows a zoom on the cubic volume containing the Milky Way progenitor halo outlined by the black square in the left image 
(46 physical kpc on a side). The grey scale background represents the dark matter projected density distribution in the Reference run. 
Overlaid circles indicate the virial radii of halos identified as being Milky Way progenitors in the Reference run (blue) and their twins in 
the Cooling (green) and Feedback (red) runs, with colours overplotted in that order (hence halos with very similar positions and radii 
in all three runs appear as red circles). A black circle is a halo in the Reference run that has no identifiable twin in the Cooling or 
Feedback run. In the right-hand image we connect the halo in the Reference run with its twin in the other runs via straight lines. The 
Milky Way progenitor halo in each run is shown in black. Most of the twins are remarkably well matched in size and position, although 
unsurprisingly, the subhalos of the Milky Way progenitor show more pronounced discrepancies between runs, especially in the central 
region of the halo. 



W e adopt a t winni n g strategy similar tolLibeskind et all 
d2010fl: IPeiranil d2010f); iRomano-Diaz et al.l l|20ld . l20Qgf iT 

ISchewtschenko fc Maccid 1 20111 ). A list of particles and as- 
sociated halo ID numbers is found for each halo in each 
simulation for each output in the Reference run. Due to the 
differences in the timestepping in each simulation there is 
typically 0.5-lMyr difference between a given Reference run 
output and a given output in the Cooling and Feedback runs. 
The list of halos in the Reference run is sorted in descending 
order of halo mass. Particles in the Reference run list are re- 
moved if they are not in halos in the Cooling or Feedback 
run. If a certain halo in the Reference run has fewer than 
50% of its particles in halos in the Cooling or Feedback run, 
it is considered to have no identifiable counterpart in the 
other simulations and is hence ignored. 

For each halo in the Reference run we calculate the frac- 
tion of its particles that belong to a halo in the Cooling and 
Feedback runs. We then select the single halo from each of 
the Cooling and Feedback runs that has the largest fraction 
of its particles in the Reference run halo and has not already 
been assigned to another Reference run halo. This provides 
a 1:1 mapping between the halos of any two simulations. 

In order to visually confirm that the twinning proce- 
dure works, we plot a map of the Reference run's projected 
dark matter density field in Fig. Q] On the same figure, we 
overplot halos in the Reference, Cooling and Feedback runs 



as circles of radius r„; r . We also link the halos in the Ref- 
erence run to their twins in the Cooling and Feedback runs 
with straight lines connecting the corresponding circles. The 
figure clearly shows that in general, the twinning procedure 
yields excellent results for most halos (the vast majority of 
circles in the left panel of Fig Q] are red). For the region 
encompassed by the Milky Way progenitor (the right-hand 
zoomed-in panel), twinning results are still quite good, ex- 
cept in the very centre where positions and sizes of sub- 
halos diverge as the non-linear nature of the system (shell- 
crossing) and the slight differences in output times between 
the runs begin to plague the comparison. 

A quantitative analysis of the twinning procedure re- 
veals that the Reference run has 96.4% of its halos twinned 
with the halos in either the Feedback or Cooling run at z=6.7 
(the final redshift for which all runs have data). If we relax 
the 1:1 mapping criterion and simply consider halos above 
the threshold where 50% of their particles in the Reference 
run also are found in halos in the other runs, 98.6% of halos 
have twins. In Fig. [2] we plot the success rate for twinning 
halos as a function of mass and redshift for the Feedback 
(solid lines) and Cooling runs (dashed lines). We find that 
for halos over 1O 9 M there is a 100% success rate for the 
twinning procedure at z=8.5 and above. At z=6.7, we find 
that the mass bins 10 8 ' 5 -10 9 ' 7 Mq exhibit a drop in twin- 
ning success rates to z=6.7 to 93% and 70% respectively. 
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Figure 2. Histograms of percentage of halos in each mass bin in 
the Reference run with twins in the Cooling and Feedback run 
against virial mass in the Reference run in Mq. The error bars 
show the sampling error on the number of un-twinned halos in 
each bin. The colours correspond to the values at different red- 
shifts (see legend). Twins in the Feedback Run are shown with a 
solid line, and twins in the Cooling Run are shown with a dashed 
line (for mass bins above 10 8 ' 6 the Cooling Run and Feedback 
run data are identical and hence the lines overlap). For halos of 
mass greater than 1O 7 M0, we find a success rate above 95% in the 
twinning procedure in most mass bins. For halos below 1O 7 M0 
(less than 200 DM particles) and at higher redshifts, this rate 
drops quite rapidly because of the lower spatial (force) resolu- 
tion in the Reference run which reduces the overall number of 
collapsed objects. For the mass bins 10 S ' 5 -10 9 ' 7 Mq , merger ac- 
tivity amongst the relatively low number of halos in these mass 
bins causes a drop in the success rate of the twinning at z=6.7 to 
93% and 70% respectively. 



This is due to the non-linearity of the N-body problem and 
merger activity as described above, as well as the relatively 
low number of halos in these mass bins. The twinning suc- 
cess rate is over 95% above 1O 7 M0 (i.e. for halos containing 
> 200 DM particles). For halos below this mass, the lower 
resolution of the Reference run causes the success rate of the 
twinning procedure to drop to between 70-90%. Note that 
only the spatial (or force) resolution in the Reference run is 
lower than in the Cooling or Feedback runs; the dark matter 
mass resolution is identical. 



2.4 Tracking high redshift galaxies down to 







The ultimate goal of this project is to compare our simu- 
lated galaxies to observed Milky Way satellites. In order to 
achieve this, we need to evolve our simulated galaxies in the 
Cooling and Feedback run to z = 0. Since it is computation- 
ally unfeasible at their nominal resolution, we instead track 
their evolution via their twin halos merger trees in the Ref- 
erence run. This determines which galaxies at high redshift 
are the progenitors of Milky Way satellites today and allows 
us to quantify how advanced satellite galaxy formation is by 
the end the epoch of reionisation. 

The fundamental assumption we make is that a halo 
which already contains stars at high redshift will still con- 
tain a galaxy at z = 0. This assumption is extremely plau- 
sible for two reasons. Firstly, even the lowest-mass halos 



are o bserved to be dark-matter dominated (|Strigari et alj 
I200&1 ). and thus we do not expect to find galaxies without 
dark matter halos. Secondly, galaxies are all predicted to be 
embedded within the inner part of the halo in which they 
form, so that the galaxy will be the last part of the halo to 
be destroyed, with ti dal stripping affect i ng the outer regions 
of the halo first (e.g. IPenarrubia et al.l (|2010l )). We further 
comment on the validity of this assumption in section 14.41 
where we identify galaxies in the Reference run at z=0 and 
locate their dark matter host halos. Finally, we also assume 
with this extrapolation technique that the dynamical fric- 
tion and tidal stripping experienced by the satellite halos in 
the Reference Run are similar to the Cooling and Feedback 
runs, i.e. that increased resolution and supernova feedback 
do not dramatically alter their efficiency. We discuss the va- 
lidity of this assumption in more detail in section 14.31 

We build the merger tree for the Reference run using the 
Branch History Method (BHM). BHM compares the sub- 
halo population of a given host halo between two snapshots, 
and attempts to optimise the tree structure to account for 
anomalies such as subhalos without an identified progeni- 
tor, or a host and subhalo switching place du ring a major 
merge r event; for details of this technique, see lTweed et al.l 
(|2009l ). As in section 12.31 we use the particle IDs to track 
dark matter particles between snapshots. For every halo in 
a given snapshot we build a list of halos in the following 
snapshot that contain particles from this halo. We then se- 
lect the halo that contains the most particles from this halo 
as its "child" halo, adopting a "one child" policy. By doing 
this, we create a halo merger tree where if halo C in output 
3 is a child of halo B in output 2 and halo B is a child of 
halo A in output 1, then halo C is also a child of halo A. 
One side-effect of this method is that if a subhalo loses more 
than 50% of its particles between two outputs, that subhalo 
is assumed to have been completely stripped by its host. To 
limit this occurrence, we use a large number of snapshots 
to build our merger tree (~ 100), so that our effective time 
resolution is roughly 150Myr. 

2.5 Resolution effects on galaxy formation 

In order to assess how reliable the Reference run is to locate 
satellite galaxies at z — 0, we first compare it to the higher 
resolution Cooling and Feedback runs in the redshift range 
where all the runs overlap. In section [231 we showed that we 
are able to very successfully match halos more massive than 
10 7 Mq between simulations. We now consider the effect that 
resolution has on star formation. 

As previously mentioned in section 12.11 the Reference 
run has the same dark matter mass resolution as the Cooling 
and Feedback runs. However, the spatial resolution, which 
determines the accuracy of both the gravitational force and 
the properties of the gas is lower; 50pc in the Reference run 
instead of 0.5pc in the Cooling and Feedback runs. The den- 
sity threshold for star formation is therefore lowered from 
10 5 atoms/cm 3 at high resolution down to 10atoms/cm 3 at 
low resolution, whilst the efficiency of star formation is pre- 
served. 

In figure [3] we compare the global star formation rate 
of halos in each of the runs. We find that before a lookback 
time of 13.1Gyr (z = 9), the Reference run's star formation 
rate is roughly half that of the Cooling and Feedback runs. 



Satellite Survival in Highly Resolved Milky Way Class Halos 7 



Redshift 

7 8 10 12 14 18 

1.00i 1 ; 1 I 1 1 1 1 




-3.00 I I 

12.9 13.0 13.1 13.2 13.3 13.4 
lookback time / Gyr 



Figure 3. Star formation rate averaged over all halos in each 
simulation. The Reference run is show in black, the Cooling run 
in green and the Feedback run in red. The vertical black line 
shows the time (z = 8.5) at which the universe is reionised in the 
simulation. Star formation in the Reference run is slightly delayed 
compared to the other runs, but catches up before reionisation. 
The jump in star formation rate at z = 9 is due to the triggering 
of a new level of refinement on the grid (from 14 to 15 levels) at 
this redshift, which allows the gas to collapse further and triggers 
star formation in all potentially star-forming halos. 
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Figure 4. Comparison of the mass-weighted stellar age for halos 
in the Reference run and their twins in the Feedback run. Halos 
that lie on the diagonal black line have the same mass-weighted 
stellar ages in both runs. Due to its lower resolution, as explained 
in section [2.51 the Reference run generally forms stars later than 
the Feedback run. However, the difference in star formation pa- 
rameters between the runs allows a modest fraction of halos to 
form stars before, especially at lower redshifts. Similar results are 
found when comparing the Reference run to the Cooling run. 

However, after 13.1Gyr, all star formation rates agree within 
30% percent. This difference of behaviour before and after 
13.1Gyr has nothing to do with reionisation, which occurs 
later on. Indeed this effect is purely numerical, and induced 
by the refinement criteria we choose to enforce. RAMSES re- 
fines the AMR grid using an octree, meaning that spatial 
resolution is a power-of-two fraction of the total box length 
(|Tevssierll2003 ). Since we specify a maximum spatial resolu- 
tion for the grid in physical parsecs, we trigger a power-of- 



two increase in resolution each time the cosmological scale 
factor has increased enough that an extra level is neces- 
sary to achieve such a resolution. In the Reference run, such 
a jum p in refinement level from 14 to 15 happens around 
z = 9. lRasera fc Tevssierl (|2006l ) demonstrate that too low a 
maximum spatial resolution delays the collapse of low mass 
haloes/galaxy disks, preventing the ISM gas density in many 
of them from crossing the star formation threshold until a 
higher level of resolution is achieved. Such a delay eventually 
vanishes when the maximum spatial resolution becomes suf- 
ficient at all times as the lack of 'step' in the star formation 
histories of the Cooling and Feedback runs on Fig [3] clearly 
shows. After a lookback time of 13.1Gyr, the star formation 
rates in all runs match well, since as discussed in section 
12.11 we select a density threshold for star formation in the 
Reference Run that best matches the star formation rates 
in the higher resolution simulations. 

We now consider the agreement between star formation 
histories of individual galaxies. In Fig. [4] we compare the 
mass weighted stellar ages of galaxies simulated at low (Ref- 
erence run) and high (Feedback run) resolution and twinned 
at 2 = 6.7. We find a pattern similar to that of the global star 
formation histories presented in Fig. [3] namely star forma- 
tion is delayed in the Reference run but converges to values 
similar to the Feedback run at a lookback time comprised 
between 13.0 and 12.9 Gyr. After this epoch there is some 
inevitable scatter due to the nonlinear nature of star for- 
mation, but this scatter is centred around the line of equal 
age in Fig. [4] The lookback time at which the ages converge 
is later than the jump in star formation due to resolution 
because the mean age is skewed by the relative paucity of 
stars formed before 13.1 Gyr in the Reference run (Fig. [3}- 
A similar result is found when comparing the Cooling run 
to the Reference run. 

It is worth pointing out that the location of star for- 
mation within a galaxy is not guaranteed to match between 
simulations. Star formation in the Cooling and Feedback 
runs is confined to regions with densities similar to molecular 
cloud cores (p > 10 5 atoms/cm 3 ), whereas in the Reference 
run star formation is allowed to occur in regions where the 
density is closer to that of typical diffuse clouds (10 at/cm 3 ). 
However, for analysing the bulk properties of satellite galax- 
ies between reionisation and z = this distinction is largely 
irrelevant. 



3 FEEDBACK IN MILKY WAY SATELLITE 
PROGENITORS 

3.1 Supernova Feedback 

We now discuss the differences between the simulations with 
and without supernova feedback, in order to better under- 
stand the role of supernovae in high-redshift dwarf galaxy 
and Milky Way satellite formation. We use the twinning 
method described in section [2~3l to match halos in the Cool- 
ing and Feedback runs. We can thus determine whether the 
net effect of including supernovae in our sub-parsec reso- 
lution simulations enhances or suppresses star formation in 
halos of various masses. We focus on halos that are captured 
by the Milky Way by z~0, and hence we can determine to 
what extent supernova feedback influences the star forma- 
tion history of Milky Way satellites observed today. Whether 
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Figure 5. Ratio of gas mass in galaxies within r v i r in the Cooling and Feedback runs, plotted against the total halo mass of their 
twin halo in the Reference run. The figure shows galaxies at the redshift of rcionization z = 8.5 on the left , and z = 6.7 on the right. 
The top plots show the ratio of total gas mass between runs, while the bottom plots show the ratio of the star-forming gas mass (i.e. 



p > lO'atoms/cnr 4 ). A black square indicates that the halo of that galaxy survives as a Milky Way satellite at 



0; grey circles 



are halos that arc completely disrupted by z = 0. We overplot as diamonds with error bars the median and interquartile range of the 
fractional differences in halo mass bins 10 7 M Q - 1O 8 M , 10 8 M Q - 10 9 M Q and 1O 9 M - 10 10 Mq . The median values lie around the 
horizontal line marking an equal ratio, with the 1O 7 M0 to 10 8 Mq bin having a ratio of ~ 0.95 and the higher mass bins having a 1:1 
ratio or higher. There is a large amount of scatter in the relative amounts of star forming gas in halos in the two simulations. It is worth 
noting that some of the halos at each redshift sampled, including all of the Milky Way satellite progenitors at z = 8.5, do not contain 
star-forming gas. This is explained as star formation occurring in bursts, with the smaller galaxies containing no star-forming gas at 
certain instants in time. 



these high-redshift galaxies become satellite galaxies of the 
Milky Way at z — or are disrupted by interactions with 
the Milky Way halo is discussed in section |4~T1 

In Figs[5]and[6] we quantify the extent to which the pos- 
itive feedback processes (metal cooling, blastwave compres- 
sion) or negative feedback processes (gas heating, outflows) 
dominate in halos of different masses. We compare the total 
gas mass, star-forming gas mass and total stellar mass in 
each halo in the Cooling and Feedback runs, using the halo 
twinning procedure described in section 12.31 Star-forming 
gas is defined as gas with a density above 10 5 at/cm 3 our 
density threshold for star formation. In each figure we over- 
plot the median and interquartile range of the fractional dif- 
ferences in the mass bins 10 7 M o - 1O 8 M , 1O 8 M - 10 9 M o 
and 1O 9 M - 1O 1O M0. For the gas masses, since there exists 
a large scatter in the results we plot the ratio for each halo 
on a log scale to highlight both large and small differences. 
For the stellar masses, since differences are smaller, we plot 
the fractional difference between the runs. In other words, 



if we denote the ratio between stellar masses by R, we plot 
the quantity (R- l)/(Q.b(R+ 1)). 

Although Fig [S] shows a large scatter in ratios of gas 
masses of twinned halos in the Cooling and Feedback runs, 
the median values lie around an equal ratio, with the 1O 7 M0 
to 1O 8 M mass bin showing ~ 5% less gas and the halos 
in the higher mass bins having a similar or slightly higher 
gas content in the Feedback run. This goes in the expected 
direction since supernovae eject gas back into the galaxy, 
causing the total gas mass to increase if this ejecta is unable 
to escape the halo. There is also considerably more scatter in 
the instantaneous star-forming gas mass results, with some 
halos in the Feedback run containing over 100 times the mass 
of star-forming gas than their Cooling run twins. This effect 
can be attributed to the enhanced metal cooling which takes 
place after the first supernovae explode in the Feedback run. 
However, on average we find a similar pattern to the total gas 
mass ratios, i.e. lower mass halos in the Feedback run contain 
very slightly less star-forming gas than their Cooling run 
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Figure 6. Fractional difference between stellar masses of galaxies in the Cooling and Feedback run (M,fb - M*Cool) / (0.5(M*fb + 
M^cool)) versus total twin halo mass in Reference Run, where M^pg is the stellar mass in the Feedback run and M^cool is the stellar 
mass in the Cooling run) . A positive value indicates that the inclusion of supernova feedback enhances star formation in the given galaxy, 
while a negative value means that supernova feedback suppresses star formation. Different panels show the values for different output 
redshifts, from z = 9 (top left) to z = 6.7 (bottom right). Black squares represent halos containing galaxies which survive as a Milky 
Way satellites at z = 0; a halo represented by a grey circle is completely disrupted by z = 0. We overplot as diamonds with error bars 
the median and interquartile range of the fractional differences in halo mass bins 10 7 M Q - 10 8 M Q , 10 8 Msolar - 10 9 M Q and 1O 9 M 
— 10 10 Mq . Star formation is slightly suppressed in low mass halos with weaker gravitational potentials. The trend is reversed for high 
mass halos. 



counterparts, and more massive halos slightly more. That 
said, at z — 8.5 the highest mass bin has a median star- 
forming gas mass that is twice as high. The large scatter 
in the amount of star- forming gas is expected; IStinson et al.l 
|2007l ) also find that star formation in their dwarf galaxies is 
quite bursty, because the instantaneous mass of star- forming 
gas can strongly fluctuate on short timescales, driven by 
catastrophic non- linear events (instabilities, mergers). 

Motivated by this result, in Fig. [6] we plot a time inte- 
grated quantity - the fractional difference between the stel- 
lar mass of each twinned halo that has formed stars in the 
Cooling and Feedback runs. Unsurprisingly, there is a maxi- 
mum 20% difference between values, much smaller than that 
for the star-forming gas mass. For similar reasons, we also 
find more scatter in stellar mass ratio of low mass halos than 
of high mass halos: the length of time that lower-mass halos 
have been forming stars is generally shorter. Therefore, they 
are more affected by temporary fluctuations in their star for- 
mation histories. As with the gas mass comparison, we find 



that the positive feedback processes outweigh the negative 
feedback processes in the highest mass bin, leading to a net 
increase in the median fractional difference in stellar mass 
of a few percent when supernova feedback is added. 



In summary, we find that the effect of feedback on the 
gas mass and star formation in a halo is complex, with lower 
mass halos being on average more affected by negative feed- 
back processes such as outflows and gas heating, and higher 
mass halos by positive feedback processes such as blastwave 
compression and metal cooling. We also find that stellar 
masses in individual twin halos can differ by up to 20%; how- 
ever, the median values only differ by a few percent in all 
mass bins, although values for low mass halos are more scat- 
tered. We therefore conclude that supernovae do not seem 
to have a significant effect on the total stellar mass of Milky 
Way satellite progenitors, at least at redshifts larger than 6. 
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Figure 7. Mean star formation histories for different halo mass bins in the Reference run (left) Cooling run (middle) and Feedback run 
(right). Mean star formation rates (SFR) in Mq /yr are given for halos of mass M < 10 8 Mq (purple astcriccs), 1O 8 M0 <M < 10 9 Mq 
(blue diamonds) and 10 9 Mq <M< 8 X 10 9 Mq (green triangles), as well as the Milky Way (orange squares) and the global SFR for the 
entire high resolution region (red solid line). The vertical line at z = 8.5 (look-back time 13.042Gyr) shows the point at which reion isation 
is turned on in the simulations. The grey region shows the detectable star formation rates as determined bv lWilkins et all d201 if) . This 
suggests that the SFR for a Milky Way-like galaxy progenitor is almost detectable at z~8 (lookback time ~13.0Gyr). We find no sudden 
drop in star formation in any mass bin after reionisation for any of the simulations. In fact, some of the star formation rates increase by 
up to 0.5 dex at reionisation. The jumps in star formation rates in the Reference run are due to the triggerin g of a new level of refinem ent 
on the grid, which allows the gas to collapse to allow star formation in all potentially star-forming regions faasera fc TevssierllioO^ . 
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Figure 8. Maximum circular velocity of satellite halos at z = versus the maximum circular velocity of their main progenitor at, from 
top left to bottom right, z = 8.5 (just prior to reionisation), z = 6.7, z = 3, z = 1. In red are the satellites that contain stars at z = 
in the Reference run; th e blue satellites remain dark. The dotted horizontal line is at 12km/s and represents the threshold given in 
lOkamoto fc Frenkl |2009) above which halos can form stars before reionisation. Like these authors, we find halos that form stars under 
this threshold to about lOkm/s; below this, no halos can form stars. Note that this threshold seems independent of redshift and that 
surviving Milky Way satellite galaxies begin to be captured by the Milky Way progenitor at z = 3. The halos with v m ax,z=o around 
80km/s arc satellites labelled 'a' and 'h' (see Fig.[9]|. 
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3.2 Reionisation Feedback 

At z = 8.5 in each simulation we include a simple instanta- 
neous, uniform heating term that represents the UV back - 
ground according to the model of lHaardt fc Madaul l|l996h . 
There are two ways such a background affects galaxy for- 
mation/evolution in our simulations: (i) it heats the ISM, 
which could prevent the gas density within galaxies from 
crossing the star formation threshold, and (ii) it heats the 
IGM, which could cut off the gas accretion onto galaxies. 
Note that we very likely overestimate these effects in the 
simulation as we neglect self-shielding wh ich is known to oc- 
cur a round densities nu < 0.1 at /cm 3 (|Susa fc Umemural 
|2004 ). On the other hand, we do not account for the local 
UV radiation (from stars within the galaxy itself) which also 
photo-ionises the ISM/IGM. Nevertheless, by looking at the 
star formation rate in halos before and after z = 8.5, we 
should be able to estimate how efficient non-local ionisation 
is at halting star formation. 

Fig. [7] shows that reionisation does not immediately 
stop star forma tion in halos a lr eady fo rming stars, in 



agreem e nt with Kitavama et al.l (|200ll); Machacek et al 



2001); iGnedin fc Kravtsov 



Okamoto fc Frenk 



20091 ); IWadepuhl fc Springell (|2010r ). Even for the lowest 
mass halos (M v i r < 1O 8 M0 ), we find that star formation 
continues after the uniform UV background is turned on. 
The free-fall time of a test particle falling from r„i r into 
one of the smallest galaxies formed at Z = 8.5 is on the 
order of 50 Myr. Hence, we conclude that star formation 
is not stopped in these halos even after ~5 free-fall times 
(the amount of time elapsed between z — 8.5 and z = 6.7). 
In other words, if reionisation does halt star formation by 
heating up the ISM or cutting off the gas accretion in halos 
that have already formed stars, it does not do so abruptly, 
but rather over a significantly extended period of time. 

Another potential effect of reionisation is to quench 
galaxy formation by preventing the collap se of gas withi n 
halos that have n ot yet formed stars (e.g. IGnedin! (|2000t ); 
ISomervillel (|2002l '); iBenson et all (|2002l M. However, the last 
z — satellite galaxy to be formed in our Reference run 
begins forming stars at z = 4.8 in a halo with M v h- = 
1.4 x 10 7 Mq , and there are 9 other satellites galaxies hosted 
by halos with a similar mass which form their first star after 
z = 8.5. Hence, whilst it is still possible that UV photoioni- 
sation has a long-term role in preventing some galaxies from 
forming, it does not seem to be able to halt galaxy formation 
entirely. 

In figure [5] we recast this statement in terms of min- 
imal circular velocity for a star forming halo, v mal , be- 
low which halos are prevented from forming stars. This al- 
lo ws us to directly compare our results to those presented 
by lOkamoto fc Frenkl <|2009h . As these authors, we cannot 
definitively conclude that this threshold arises entirely be- 
cause of reionisation or the general inability of halos below a 
v mal of lOkm/s to cool and form stars by z = 0, since we do 
not run a simulation without reionisation. However, we note 
that reioinisation in our simulation occ urs instantaneously 
at z = 8.5 (close to the value of z — 9 of lOkamoto fc Frenkl 
|2009h ) and that, in stark contrast the threshold of v ma x ~ 
10 km/s seems independent of redshift. Indeed, it remains 
quite constant both before and after reionisation has oc- 
curred, which leads us to argue that reionisation cannot play 



an important role in setting its level and only sustains it, in 
the best of cases. 



4 MILKY WAY SATELLITES TODAY 
4.1 Tracking satellite galaxies down to z = 



In section [2741 we discussed the techniques used to track the 
galaxies formed in the Cooling and Feedback run down to 
z = using the Reference run. We now analyse the results 
of this tracking. We resolve about 6630 such halos. Of these 
6630, 394 survive as subhalos of the Milky Way halo at z = 
0. In table[2] we list the number of galaxies formed by various 
redshifts from z = lltoz = 6.7in the Cooling and Feedback 
runs that survive to z = 0. It is apparent from the table 
that by this redshift, we have not formed enough satellite 
galaxies to m atch even the population of pre-SDSS satellites 
(Mated Il998h . However, it is also clear that satellite galaxy 
formation continues after reionisation; three galaxies that 
end up as satellites of the Milky Way at z — are formed 
between z = 8 and z — 6.7. 

In fact, we find that the Milky Way satellite halo that 
began forming stars latest in the Reference run does so at 
z = 4.8, or a lookback time of 12.4 Gyr, after reionisation is 
complete. This is illustrated on Fig. 1111 where we plot the 
stellar mass of each Milky Way satellite galaxy in the Refer- 
ence run at z = against the age of their oldest star particle. 
As we note in section 12.21 we are unable to identify galax- 
ies with a stellar mass below 3.5xlO 5 M0 in the Reference 
run. Hence we cannot discount the possibility that the Cool- 
ing and Feedback runs might form more galaxies with lower 
masses that survive as Milky Way satellite galaxies. All we 
can conclude is that every galaxy above this mass thresh- 
old that survives as a Milky Way satellite in the Cooling 
and Feedback runs (through the twinning procedure) also 
survives as a satellite in the Reference run. However, it is 
possible that the trend of lower mass satellite galaxies form- 
ing at lower redshifts continues in the Cooling and Feedback 
runs, where stellar masses is better resolved. 

Key to the survival process of satellite galaxies is the 
mass stripping they undergo as a function of time. We visu- 
alise this in Fig [9] where we identify which halos containing 
stars at z — 6.7 survive to become Milky Way satellites at 
z = and follow their dark matter particles through cosmic 
time. We locate these particles in outputs of the simulation 
at z = 3, 1 and 0, and overplot them on top of their respec- 
tive underlying density fields, colour coding them according 
to the redshift at which the halo they belong to is captured 
by that of the Milky Way. We find that halos captured be- 
fore z = 1 experience significant disruption despite the core 
of the halo surviving (halos 'c', 'e', 'f and 'g' in Fig. [9}. (By 
contrast, we find that the only stars formed by z = 6.7 that 
are stripped from their galaxies are from satellites captured 
before z = 3). The surviving satellite that is captured just 
after z — 3 (halo 'g') arrives in the Milky Way halo as a sub- 
halo of another halo, which is subsequently disrupted by the 
Milky Way. Hence this halo has already experienced some 
stripping by z — 3, as shown in Fig. [5] 

In Fig. 1101 we investigate satellite survival to z = in 
more detail. We compare the redshift at which a halo is cap- 
tured by the Milky Way halo against its mass at capture in 
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Figure 11. Time of first star formation for satellite galaxies that 
survive to z = (including the MW itself) in the Reference run 
against their stellar mass at z = 0. The region shaded in grey 
represents the lookback times which have not been simulated in 
the high resolution Cooling and Feedback runs. This illustrates 
that satellite galaxy formation is incomplete until z = 4.8 (i.e. 
after a lookback time of 12.4 Gyr). Note that this plot only shows 
galaxies in the Reference run. 



two simulations: the Reference run and the matching Adi- 
abatic run. The only difference between these two runs is 
that an extra right hand side 'sink' term is included in the 
energy equation of the gas in the Reference run to model 
losses due to radiative cooling, as well as star formation (see 
section [2HJ. We find that out of all the halos which survive 
to z — in these two runs (6 in the Reference run, 20 in the 
Adiabatic run), only one is captured at z > 3 in the Adia- 
batic run and none in the Reference run. Moreover, higher 
mass halos (M v i r > 10 9 ~Mq) only survive if they are cap- 
tured later so that the highest mass satellites at z = are 
systematically the ones that are captured last. It thus ap- 
pears that the inclusion of radiative cooling in simulations 
of Milky Way-like galaxies has a dramatic impact on the 
survival of the satellite galaxies. We discuss the reasons for 
this discrepancy in section 1431 

their pure dark matter counterparts. This is explained 
by two effects: (i) it is more difficult to strip mass from 
satellites as their central density increases and (ii) the cus- 
piness (and central density) of the host halo is increased. It 
is interesting to note that effect (i) could in principle lead 
to the opposite effect, i.e. an increase in the lifetime of the 
satellites as it makes them more concentrated and thus more 
resistant to tidal disruption, but we find that reduction in 
dynamical friction time scales due to mass increase domi- 
nates. Of course these conclusions could be altered if a sub- 
stantial mass of baryons was ejected out of the satellites, 
to the poin t where the trend that we m easure could even 
be reversed. IPontzen fc Governatol (|2012l ) argue that such a 
reversal is plausible, though it would require a significantly 
more efficient feedback mechanism than the one we observe 
in section [3~T1 



4.2 Dark Matter Satellite Halos 

In this section, we consider the population of dark matter 
subhalos of the Milky Way in our runs at z = 0, comparing 
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Figure 12. Cumulative maximum circular velocity (Vmax) func- 
tions comparing our results to high resolution dark matter only 
simulations of MW-like objects at z = published in the 
literature. The Reference run data is shown as a thick solid 
line, and the Dark Ma tter run as a thick dashed curve. Via 
Lactea II (VLII) data llDiemand et alj|2008h is overplotted as 
a thi n dotted curve. A fit to the Aquarius data l lSpringel et al.l 
l20Qgf) is plotted as a t hin dashed l i ne. T he empirical formula 
for N(<Vmax) given by iReed et al. | J2005I) is p l otted as a thin 

' Mated Jl998V) : iBekki fc Chibal 



solid line. Observational data from 



ll2005l) : lBekki fc Stanimirovid {20091) : I Wolf et aU ||2010|) is show n 
as diamonds. In each case, Vmax is given by max(- v / GAf (< r)/r), 
where M is the mass inside rso, the radius at which the density 
exceeds 50p cr ;t. 



and contrasting their properties with similar dark matter 
simulations which exist in the literature. 

For this purpose, we run a fourth simulation, the Dark 
Matter run, which is a pure dark matter version of the Ref- 
erence run (see section I2.1|l . The results of this simulation 
are described in more detail in section [4731 W e use it to com- 
pare our results directly with the Aquarius dSpringel et all 
l2008h and Via Lactea II (|Diemand et al.ll2008h simulations, 
which are the most resolved dark matter only simulations 
of MW-like objects available to date. The most basic com- 
parison, a cumulative maximum circular velocity function, 
is presented in Fig. 1121 In this fig ure, we also overp lot the 
empirical prescription proposed bv lReed et al.l (|2005l ) for the 
Milky Way halo in our Dark Matter run at z=0. We find that 
our Dark Matter run satellite halo data is well represented 
by this prescription, but that the Reference run predicts sig- 
nificantly fewer satellites at the low velocity end (between a 
factor 2 and 3 for Vmax < 20km/s), and more massive satel- 
lites {Vmax > 30km/s). We discuss the impact of simulation 
physics on the maximum circular velocity function in the 
next section (|4.3[) . 

From figure [121 it is apparent that, whilst the cumula- 
tive maximum circular velocity function of our Dark Mat- 
ter run has the same shape {N(> Vmax) oc Vmax) as that 
measured in both the Aquarius and Via Lactea II simu- 
lations, its normalisation is more than an order of magni- 
tude lower. This large discrepancy can almost entirely be 
attributed to our choice for th e mass of the MW host halo 
since our agreement with the IReed et all (|2005l ) prescrip- 
tion is quite reasonable (better than 20 %). Our host halo 
has a Vmax of 126km/s in the Dark Matter run as op- 
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captured by MW 

Figure 9. Tracking the location of dark matter particles in MW satellite progenitor halos from 2 = 6.7 to z = in the Reference run. 
Each image is a projection of a cubic volume of length shown on the x-axis in physical kpc. Blue circles represent the virial radius of each 
halo tracked; the MW halo is shown as a black circle. The image at z = lies largely inside the Milky Way virial radius. The colour of 
the particles in a halo represents the redshift at which the halo is captured and becomes a subhalo of the MW (see the colour bar). Halos 
captured before 2 = 1 exhibit significantly more stripping at z = than halos captured after 2 = 1. Halo 'g' (in orange) is captured and 
partially stripped by another halo which, in turn, is captured and completely disrupted by the Milky Way between 2^3 and 2 = 1, 
while halo 'g' itself survives as a luminous satellite galaxy of the MW at z = 0. 
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Figure 10. Survival of halos as MW satellites at z = 0. Both panels show the redshift at which a halo becomes a MW sub- halo, 2 cap ture, 
against the sub-halo mass when this capture happens, M cap t U re; on the left is the Reference run, and on the right is the Adiabatic run. 
Black squares survive as MW satellite halos at z = 0. Note that they can be completely disrupted as long as they become part of another 
MW satellite halo by z = in the disruption process. Grey circles are completely disrupted and simply become part of the diffuse MW 
halo by z = 0. Halos captured by the MW host before z = 3 do not survive to z = 0, with one exception in the Adiabatic run. Less than 
half of the low-mass halos (M v i r ^ 10 9 Mq) captured after z = 3 are able to survive to z = in the Reference run, whereas the vast 
majority of them survive in the Adiabatic run. 
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Table 2. Table of the fate of galaxies formed between 2 = 11 and z = 6.7. The 6 columns are, from left to right: (1) the redshift at 
which the stellar population is sampled in the Reference, Cooling and Feedback runs ('high z'); (2) the simulation name; (3) the number 
of galaxies which become satellites of the MW between the sampled redshift and z = not including the main MW progenitor halo (in 
brackets, the number of those galaxies whose twin halos in the Reference run also contain at least a galaxy); (4) the number of galaxies 
surviving as MW satellites at z = 0; (5) the satellite progenitors destroyed by mergers with other satellite progenitors - the figures in 
brackets show the number of objects taking part in mergers at high z, followed by the resulting number of objects after the satellite 
progenitor-satellite progenitor mergers at z = 0; (6) the number of galaxies merged with the MW and destroyed between high z and 
2 = 0. See section [2.41 for a description of how these numbers are calculated. We find that the large majority of the halos containing 
galaxies captured by the MW by z = merge with it and are destroyed, with two galaxies merging with each other before being captured 
by the MW and becoming a satellite galaxy. More mergers of MW progenitor galaxies are found, but these are all completely disrupted 
and destroyed after capture by the MW. We find that satellite galaxy formation in the Reference run is not complete by the lowest 
redshift reached by the high resolution runs (z=6.7). By looking at the ages of the star particles in satellite galaxies at z = in the 
Reference run, we find that satellite galaxy formation continues until at least z = 4.8 (see Fig. Illll . 
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posed to ~200km/s in the Aquarius or VLII simulations, 
i.e. 40% less than either of these simulations. Note that in 
figure 1121 we plot all subhalos inside rso (denoting the ra- 
dius at which the density is above 50p cr it) as opposed to 
T200, which we u se in the res t of th is paper. This is to allow 
comparison wi t h iReed et all (|2005l ); IPiemand et all (|2008h ; 
ISpringel et al.l (|2008l ). By contrast, the Reference run lies 
below t he empirical formula, which we comment on in sec- 
tion 14.31 ISpringel et al.l (|2008l ) r eport that their si mulations 
overshoot the fitting formula of IReed et al.l (|2005l ) by a fac- 
tor ~ 3 which they argue most likely arises from a systematic 
effect in the numer ical technique used t o perform the runs. 
On the other hand, iMadau et al.l (|2008l ) suggest it is due to 
the different (WMAP 1 instead of WMAP 3/5/7, i.e. dif- 
ferent normalisation, erg, and/or tilt, n s , of the power spec- 
trum) cosmology the Aquarius simulations employ. In any 
case, the remarkable conclusion that we draw from this com- 
parison exercise, is that for a MW host halo with V ma a: = 126 
km/s in the Dark Matter run as opposed to ~ 210 km/s in 
the Aquarius or VLII simulations, i.e. a difference of about 
65 %, one gets a suppression in the number of satellites by 
about a factor of 10, which is enough to match the observed 
abundance of MW satellites with Vmax between 10 and 30 
km/s. Now, there is still an ongoing debate as to what the 
exact mass of the Milky Way halo i s (e.g . Battaglia et a! 



2005 1 ) ; iKarachentsev fc Kashibadzej (|2005l ); IWatkins et al 



20101) ). Our simulated halo admittedly lies at the very low 
end of the estimated range of values (4.32 x 10 n MQ within 
195 kpc), and Aquarius and Via Lactea II somewhat on the 
high side (1.85 xlO 12 M within 245 kpc). It is not the 
purpose of the present paper to constrain the Milky Way 
mass, but only to illustrate how the uncertainty in the mass 
of the Milky Way halo and the inclusion of baryonic physics 
translate into an uncertainty in the number and properties of 
MW satellites one predicts. Thus we simply remark that the 
observed abundance of these satellites, taken at face value, 
seems to favour a less massive MW halo. 



4.3 The Effect of Baryonic Physics on Satellite 
Galaxy Survival 

As described in section 12.11 we run three simulations to 
z = 0, the Dark Matter Run, the Adiabatic Run and the 
Reference Run. We also follow galaxies formed by z = 6.7 
in the Cooling and Reference run to z = by using the Ref- 
erence run's merger tree. Only the Reference Run includes 
star formation up to z~0. The first thing to note is that 
the Dark Matter run produces more (sub-)halos (~ 20% 
more) than any of the runs containing baryons; we present 
exact numbers in table [3] These results are found both at 
z=7 and z=0. It thus appears that pure dark matter simu- 
lations, which do not include baryonic pressure forces, allow 
for the more efficient collapse of halos than simulations that 
contain baryons. 

Secondly, as described earlier in section 14.11 we find 
fewer subhalos in the Reference run than in the two other 
runs. In figure 1141 we quantify this effect by plotting the 
number of high redshift galaxies surviving to become satel- 
lite galaxies of the Milky Way at z — (see sections 12.41 
and 14.11 for details on the tracking of satellites from high 
to low redshift and table [2] for numbers specific to the Ref- 
erence run). We find that the Reference run contains far 
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Figure 13. Density profiles of the main Milky Way-like halo up 
to r v i r at z = for the Reference (thick solid line), Adiabatic 
(dotted line) and Dark Matter (dashed line) runs. We plot total 
density, i.e. dark matter, gas and stars in the Reference run and 
dark matter and gas in the Adiabatic run; the dot-dashed line 
represents the DM density in the Reference run. We overplot a 
vertical dotted line at 200pc, equivalent to AMR level 16 at z = 0, 
i.e. 2 levels below the highest resolution reached by the runs to 
indicate the scale below which the gravitational force is underes- 
timated, and lines at r2oo, used in the bulk of the text and rso, 
used for comparison purposes with figure 1121 Power law profiles 
scaling like r _1 and r~ 2 and a NFW profile fit to the DM run 
(thin solid line) are also overplotted. 
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Figure 14. Number of satellite galaxies at z = formed at or 
before a given redshift. Satellite galaxies are defined as sub-halos 
whose progenitors at high redshift have a twin in the Feedback 
run that contains one or more galaxies (see section 14.111 . Each 
line represents a different run used to track halos down to z = 0. 
The solid line shows the Reference run, the dashed line shows 
the Dark Matter run and the dotted line shows the Adiabatic 
Run. The inclusion of gas cooling reduces the number of galaxies 
at high redshift that survive down to z = by a factor ~ 3 to 
10 depending on redshift. As discussed in the text, the central 
density profile of the Milky Way at z = plays a key role in 
governing the survival of satellite galaxies. 
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Redshift Halo Type Dark Matter Adiabatic Reference Cooling Feedback 



Independent Halos 2084 1706 1731 

Subhalos 2730 2418 2209 



Independent Halos 2405 1997 1997 2142 2104 

Subhalos 598 449 343 445 390 

Table 3. Total number of independent halos and subhalos (dark and luminous) in the Dark Matter, Adiabatic and Reference runs at 
two sampled redshifts. Runs containing baryons compare favourably, whereas the Dark Matter run contains significantly more halos for 
a given redshift. There are also fewer subhalos in the Reference run than either of the other two runs, as explained in section 14.31 



fewer (factor 3 to 7 depending on redshift) surviving satel- 
lite galaxy sub-halo hosts than the Dark Matter or Adiabatic 
runs. Similarly, figure ITTJ1 shows that while in the Adiabatic 
run, all halos captured by the Milky Way after z — 3 sur- 
vive to z = 0, in the Reference run only the less massive 
halos survive to z = 0. Finally, we note that in figure 1121 
the Reference run subhalo cumulative maximum circular ve- 
locity function at z = lies below the Dark Matter run's. 
There is hence a body of evidence to suggest that including 
gas cooling in our simulations of a Milky Way-like halo has 
significantly altered the population of halos, both dark and 
luminous, surviving to z — 0. 

This discrepancy can be explained as the result of in- 
creased dynamical friction on the more massive infalling 
satellites, which is exacerbated by differences between the 
simulations in the density profile of the Milky Way. Indeed, 
the Chandrasekhar formula, which correctly encapsulates 
the basic physics of dynamical friction, states that the fric- 
tion force is proportional to the mass o f that object and on 
the d ensity of the surrounding medium (|Binnev fc Tremainel 
120081 ). Figure [13] shows that the total (gas & DM) density 
profiles of the main MW halo at z = are very similar be- 
tween the DM and Adiabatic runs, with the Adiabatic run 
having a slightly lower density overall. However, the profiles 
diverge at around lOkpc; at a radius of 2kpc, the density in 
the Reference run is several times higher than in the other 
runs. This is due to the fact that gas in the Reference run is 
able to cool, and hence the density of the halo within a ra- 
dius of lOkpc is significantly higher in this run, since cooled 
baryons ar e able to condense at the c entre, pulling DM along 
with them lFlores fc Primackl (|l994l ). 

We analyse the halos containing satellite galaxies that 
survive to z = in the Adiabatic run but not the Refer- 
ence run and find the following. All of these halos pass close 
to the Milky Way's centre in both the Reference and Adia- 
batic runs. In addition, when close to the centre, the halos' 
velocities are highly radial, suggesting that these subhalos 
have experienced sufficient dynamical friction and lose most 
of their angular momentum. Depending on their mass and 
density, the satellites survive 1 to 7 passes before their orbit 
decays to the point where they merge completely with the 
Milky Way halo in the Reference run. In the Adiabatic run, 
where the central density of the Milky Way is much lower, 
the orbits decay much more slowly, and as a result these 
objects are not destroyed by z = 0. By comparison, the sur- 
viving satellites in the Reference run do not pass close to 
the centre of the Milky Way, and thus survive to z = 0. 

To conclude, even if our sub-halos had exactly the 
same mass in the DM and Reference runs, a higher cen- 
tral density of the host halo would cause their orbits to 
decay faster, thus decreasing their survival time, as found 



bvlRomano-Dfaz et al.l (|2010l ) and lSchewtschenko fc Macciol 
120111 ). This effect is exacerbated as the satellites in the 
Reference run are also more concentrated, rendering them 
less prone to tidal disruption at large radii, retaining more 
mass on their way to the centre of the host. This higher 
mass makes them more susceptible to dynamical friction. It 
should be noted that we use an AMR code, as opposed to 
SPH as previous authors do. One consequence of this is that 
gravitational force resolution is dependent on the resolution 
of the grid structure. However, due to our relatively high res- 
olution (50pc in the Reference run), our force resolution in 
the centres of the MW halo and its satellites is competitive 
with, if not better than, that of previous authors. We thus 
believe that the main caveat of our Reference run is that it 
does not include a feedback model, although we point out 
that standard supernovae feedback, as modelled in our high- 
resolution Feedback run, is extreme ly unlikely to reverse 
the si tuation significantly. However, iPontzen fc Governatol 
|2012l ) have recently argued that the injection of energy from 
supernovae in the centre of dwarf galaxies at 4 > z > 2 can 
dramatically alter their halo density profiles. Whilst we are 
not able to confirm this effect with our own set of simula- 
tions, should it prove to be able to lower the central density 
of the Milky Way halo as well, we predict that more satel- 
lites will survive to z = 0. However, whether this number 
will match that measured in pure DM simulations or will 
still reflect a significant suppression of satellites is likely to 
depend on the details of the numerical implementation of 
the feedback processes. 



4.4 Properties of Satellite Galaxies at z = 

In this section we compare the properties of galaxies found 
at z — in our simulations to observations of known satel- 
lite galaxies of the Milky Way and M31. To do this, we 
make use of the Reference run results measured at z = 
directly rather than relying on the satellite tracking algo- 
rithm described in section 14.11 In this run, we find fifteen 
satellite galaxies using a direct galaxy identification method 
described in section [2~2l Four of these (labelled i, I, r and q 
on Fig. I16p do not have associated dark matter host halos 
identified by our halo detection algorithm. Three of these, 
i, I and q, once had hosts that fell below our halo mass 
resolution threshold of 2.2 x 10 6 M by z=0. The host of 
the fourth, r, merged with another halo at z=5.3, which then 
merged with the Milky Way at z=2.3. We would thus require 
a higher dark matter resolution to comment on whether 
satellite galaxies are capable of losing their dark matter host 
halos without themselves being destroyed by stripping. We 
further note that including a DM sub-halo host with a mass 
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Figure 15. Star formati on histories of satelli te gal axies in the Refere nce run (black squares) compared to values deduced from observations 
of Milky Way satellites iOrban et al. | feOOSl) and lWolf etafl j2010l) : grey circles). Each panel shows the fraction of stellar mass formed 
between z = and the lookback time indicated on ordinate axis against the half stellar mass of the galaxy at z = 0. Stellar fractions 
and masses in each plot match well with the values estimated from observations, suggesting that the star formation histories of satellite 
galaxies in the Reference run do not wildly differ from those of observed Milky Way satellites. 



just below our mass resolution threshold (2.2 x 10 6 Mq) 
would not affect the results we present here. 

Figure [16] shows a plot comparing the half stellar mass 
against half-light radius of our simulated satellite galax- 
ies (grey sq uares), to t he sat ellite galaxies of the Milky 
Way (from IWolf et al.l (fold)) ( white circles) a n d M3 1 
(from ICollins et all J201ll. l201Ch: iFerguson et all (|200d ): 
Irwin et all j200Sl): llrwin fc Collins] (l201ll, priv. comm.); 
Kaliraietal] i201Ch: iLetarte et al. J2009I); iMartin et all 
J2009I); iMcConnachie et alj J2008I); iMorrison et all {200$ ); 
Pustilnik et all (|2008l ); iRichardson et al.l (|201ll ) ) (white 
squares). The population of simulated galaxies lie in the 
same region as the observed satellites except for the two 
galaxies labelled a and h on the figure, which lie well above 
the observational data. These are the most massive galaxies 
in our sample, which have recently been captured by the 
Milky Way host halo and thus have not experienced signif- 
icant stripping by z — 0. We note that beyond their larger 
stellar half ma ss, which c a n be somewhat overlooked since 
the data from IWolf et all (|20ld ) do not contain the LMC 
or SMC, these simulated galaxies seem too compact when 
compared to the rather well defined observational relation 
linking satellites size and mass. While resolution undoubt- 
edly plays an important role in getting accurate estimates of 



the sizes of simulated objects, it is nevertheless striking that 
satellite galaxies which lie on top of the observed results all 
exhibit large tidal tails (figure [9]) and thus evidence of tidal 
stripping, compared to the two most massive galaxies ('a' 
and 'h' in figure [5]), which enter the Milky Way halo much 
later, after z=l. This is relevant because the LMC and SMC 
are believed to have e ntered the Milky Way later than the 
other satellite galaxies (|Besla et al.ll2O10l ). which can explain 
the morphological differences between them. 

Moreover, when we compare star formation histories of 
satellite galaxies in the Reference run with those derived 
from the analysis of colour-magnitude diagr ams (CMDs) 
of Mi lky Way satellite galaxies observed by lOrban et al.l 
(|200Sh we find reasonable agreement. This is demonstrated 
in Fig. 1151 where we plot the fraction of stars formed af- 
ter lookback times of lGyr, 2Gyr, 5Gyr and lOGyr both in 
simulated (grey squares) and observed (white circles) galax- 
ies. It is somewhat reassuring that the star formation his- 
tories of our simulated satellite galaxies are broadly correct 
since we have previously shown (Fig [4]) that they are fairly 
independent of resolution. Perhaps more importantly, this 
agreement also suggests that feedback, whatever its form 
and origin, cannot drastically alter the star formation histo- 
ries of these galaxies: models where significant feedback at 
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Figure 16. Properties of satellite galaxies surviving down to 
2 = 0. The top panel shows half stellar mass against half- 
stellar-mass radius for satellite galaxies in the Reference run at 
z = (large black squares), compared to half-light radii cal- 
culated for satellites of the Milky Way and M31 (small grey 
squares and small grey circles respecti vely). Data for th e satel- 
lites of the Milky Way are taken from IWolf et al, I tepid) . Data 



l2010t);lFereuson et al. 


J2000l);llrwin et al. 


d2008t); Irwin & Collins! 


ll201ll. priv. comm.); 


Kalirai et al.l (12010 


); iLetarte et al.l J2009I); 


iMartin et alj J2009I): 


McConnachie et al. 


(2008 


); Morrison et al.l 


(2003); Pustilnik et al. (2008); Richardson et al 


J201lh . The bot- 



torn panel shows the spatial location of the satellite galaxies and 
their host DM sub-halos inside the Milky Way virial radius (red 
circles indicate galaxy radii, while blue circles represent sub-halo 
radii). The black circle shows the Milky Way virial radius. The 
background image is a cubic projection of the stellar density field. 
The two satellite galaxies (labelled a and h) which lie above the 
observed data in the top panel are galaxies that have been lit- 
tle affected by stripping due to their late capture by the Milky 
Way host. Note that the Milky Way data does not include the 
Magellanic clouds or ultrafaint satellites. 



z > 1 completely quenches late star formation are clearly 
ruled out by the observational data. This makes it all the 
more challenging for stellar feedback to soften cusps of dark 
matter sub-halos and certainly favours a rapid, irreversible 
m echanism, very localised in tim e such as the one suggested 
bv lPontzen fc Governatol (|2012T l. 



5 DISCUSSION AND CONCLUSIONS 

The work discussed in this paper has made use of the Ssr^sS 
suite of high (a few tens of parsec) to ultra-high (sub-parsec) 
resolution cosmological hydrodynamic re-simulations to in- 
vestigate the effect of baryonic physics on the evolution of a 
'Milky Way' and its satellite galaxies. Whilst various other 
authors have simulated satellites of Milky Way-like galaxies 
down to z = 0, ours are the first to reach sub-parsec res- 
olution down to the end of the epoch of reionisation. The 
motivation for this was to analyse in detail the evolution of 
satellite galaxies around the epoch of reionisation, which has 
been posited as a mechanism for suppressing star and galaxy 
formation in dwarf halos and hence for shaping the popula- 
tion of satellite galaxies that we observe around the Milky 
Way and M31. To the best of our knowledge, we are also 
the first authors to use an adaptive mesh refinement (AMR) 
technique to study the evolution of these satellite galaxies 
down to 2 = (other studies thus far either used smoothed 
particle hydrodynamics (SPH), or ended their simulations 

at higher redshifts). 

iGuo et al.l 



In agreeme nt wit h e.g. 



2010) 



2010), 



iRicotti fc Gnedinl (|2005l) ; IWadepuhl fc Springell 
we find that reionisation appears not to efficiently stop star, 
or even galaxy formation. Instead, we find that satellite 
galaxy formation continues down to at least z = 4.8 in 
our lowest (~50 pc) resolution simulation. The number 
of luminous satellite galaxies formed before reionisation 
(z = 8.5 in our case) is found to be far lower than the 
number of observed Milky Way satellites. The s e results are 
consistent with e.g. lOkamoto fc Frenkl (120091 ); iHoeft et al.l 
(j2006T) : like these authors, we find that, down to the end 
of the reionisation era, there exists a threshold in v max 
of about 10 km/s below which halos remain dark, never 
forming stars. This threshold persists at later redshifts, 
i.e. well after reionisation has ended. This is probably due 
to the fact that efficient atomic gas cooling increases halo 
central densities and hence Vmax, separating halos that 
can cool gas and form stars from those that cannot. There 
are, however, at least two major limitations in our work 
that prevent us from commenting further. First, we do not 
run a simulation with self-consistent star formation and 
ionisation, and hence we cannot quantify the precise effect 
of UV photoionisation on star formation in galaxies already 
forming stars. Secondly, our model of reionisation consists 
of a uniform background which neglects the effect of gas 
self-shielding from external photoionisation sources. Future 
studies which include ionising photon radiative transfer 
and hence self-consistent re-ionisation, should allow us to 
determine whether proximity effects and/or self-shielding 
significantly alter our conclusions, but this does not seem 
likely. 

The effect of supernova feedback on the gas and stel- 
lar mass of high redshift (both pre- and post-reionisation) 
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galaxies seems to be quite stochastic. Indeed, our results ex- 
hibit a large scatter, even though one might argue we detect 
a slight systematic trend of star formation being enhanced 
by feedback as galaxy mass increases. A key feature of our 
model of supernova feedback is that it consists of star parti- 
cles (with a standard Salpeter IMF) injecting mass, energy, 
momentum and metals into the surrounding gas according 
to a Sedov blast wave solution deposited onto the grid lOMyr 
after they have formed. Together with our sub-parsec reso- 
lution, this means that we should be able to track the effect 
of supernovae explosions fairly realistically on scales typi- 
cal of small molecular clouds. Outflows and heating from 
these supernovae reduce the amount of gas available for star 
formation, but blast wave compression and an excess of ra- 
diative cooling due to the injection of metals into the ISM 
can potentially increase it. Our results show that, on av- 
erage, adding supernovae feedback does reduce the gas and 
stellar mass of galaxies in halos below 10 9 Mq (negative feed- 
back) but increases it for galaxies hosted by halos above 
10 9 Mq , where the deeper potential and extra metal injec- 
tion negates the impact of outflows (positive feedback). In 
any case, the stellar masses of individual galaxies are only 
changed by maximum 10-20 percent either way by super- 
nova feedback, and one has to invoke a much larger ener- 
getic input (for example a top-heavy IMF, or a large fraction 
of hypernovae) and/or one that is impervious to radiative 
losses (e.g. some kind of 'turbulent' energy) to overturn the 
situation in favour of negative feedback. 

When we use lower resolution (50 pc) runs to track 
the descendants of galaxies in the high resolution (0.5 pc) 
Cooling and Feedback runs to z = 0, we find that very few 
of the galaxies which are captured by the Milky Way pro- 
genitor halo survive to z = 0. In fact, independent of the 
run used to perform the tracking (PureDM, Adiabatic or 
Reference), no satellite galaxy captured by the Milky Way 
progenitor before z = 3 survives as Milky Way satellite at 
z = 0. However, the main difference is that all galaxies (sub- 
halos) captured after z = 3 in the Pure DM or Adiabatic 
runs survive, while only a small fraction of these does so in 
the Reference run. This is caused by the much higher cen- 
tral density of the Milky Way halo (and sub-halos) in the 
Reference run (the only one to host a MW galaxy), which 
significantly shortens the dynamical friction timescales for 
satellites t o spiral to the centre o f the halo and experience 
disruption. iLibeskind et al.l (|2010T i suggest that satellites in 
simulations including gas cooling (rather than pure dark 
matter) experience lower mass l oss, although they are more 
radially concentrated. However, iRomancn Di'az et all (|201Ch 
find, as we do, that inclu ding baryon physics redu c es the 
survival time of satellites. ISchewtschenko fc Macciol (12011 



explai ns this discrepancy by noting that ILibeskind et al.l 
\201(t ) do not measure satellite survival but rather mass 
loss over time, whereas Romano- Diaz et al.l |2010l ) and 
ISchewtschenko fc Macciol ( 20111 ) find that it is in the cen- 
tre of the halo that the satellites experience the most mass 
loss, and hence the centre of the host is where the survival 
of satellites is determined. 

We believe that the work presented here, which makes 
use of a completely different simulation technique and im- 
proves on the resolution of these previous studies, sheds a 
useful light on the issue. In particular, we emphasize that 
when we compare properties of the remaining satellite galax- 



ies at z = in the Reference run to their observational 
equivalents for the Milky Way or M31, we find that the star 
formation histories, stellar masses and radii are in reasonable 
agreement. Only the more massive, recently captured satel- 
lites are found to be too compact when compared to Milky 
Way and M31 dwarf spheroidal satellites. This puts rather 
tight constraints on the feedback mechanisms (amount of en- 
ergy, duration, timing) required to soften the cusps of dark 
matter halos as they cannot have a major impact on these 
properties. 

Finally, beyond the importance of the role played by 
baryonic physics in determining the number of satellites of 
MW class halos, it is worth noting the extreme sensitivity of 
this number to the circular velocity of the host halo. Whilst 
all simulated CDM halos in the lite rature match the shape 
of the iV sa t cx Vmam relation from iReed et al.l (|2005T l. not 
all of them agree on the constant of normalisation. Indeed, 
our Pure DM run agrees with the normalisation given by 
IReed et al.l (l2005h at the 10-20 % level whereas the Aquar- 
ius ( Springel et alj|2008h simulations overshoot it by a fac- 



tor '^_3 J _and_by about 30 % more than the Via Lactea 
II (|Diemand et alj|2008l ) simulation. Curiously enough, the 
cause of this significant discrepancy is still unresolved, even 
though the d ebate has received a lo t of attention lately in 
>aper s s uch as Vera-Ciro et al.l (|201lh ; feovlan-Kolchin et al.l 
2012t ); IWang et al.l (|2012T ). Irrespective of this disagree- 
ment, our results, especially when baryonic physics is taken 
into account, argue in favour of a less massive (5 xlO 11 
Mq < M vlr < 10 12 M ) MW halo. Note that if the like- 
lihood for a satellite galaxy to survive down to z = is 
suppressed by even half as much as we find in this work 
when feedback is properly incorporated, models relying on 
N-body dark matter-only simulations, such as semi-analytic 
models (SAMs) should be revised accordingly. 
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